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ABSTRACT 

Recently, we have shown that if the ISM is governed by super-sonic turbulent flows, the excursion-set 
formalism can be used to calculate the statistics of self-gravitating objects over a wide range of scales. On 
the largest self-gravitating scales ("first crossing"), these correspond to GMCs, and on the smallest non- 
fragmenting self-gravitating scales ("last crossing"), to protostellar cores. Here, we extend this formalism 
to rigorously calculate the auto and cross-correlation functions of cores (and by extension, young stars) as 
a function of spatial separation and mass, in analogy to the cosmological calculation of halo clustering. We 
show that this generically predicts that star formation is very strongly clustered on small scales: stars form 
in clusters, themselves inside GMCs. Outside the binary-star regime, the projected correlation function 
declines as a weak power-law, until a characteristic scale which corresponds to the characteristic mass 
scale of GMCs. On much larger scales the clustering declines such that star formation is not strongly 
biased on galactic scales, relative to the actual dense gas distribution. The precise correlation function 
shape depends on properties of the turbulent spectrum, but its qualitative behavior is quite general. The 
predictions agree well with observations of young star and core autocorrelation functions over ~ 4 dex in 
radius. Clustered star formation is a generic consequence of supersonic turbulence if most of the power 
in the velocity field, hence the contribution to density fluctuations, comes from large scales ~ h. The 
distribution of self-gravitating masses near the sonic length is then imprinted by fluctuations on larger 
scales. We similarly show that the fraction of stars formed in "isolated" modes should be small, < 10%. 
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1 INTRODUCTION 

A central tenet of our understanding of star formation is that most 
stars form in clusters (e.g. Lada & Lada 2003; Portegies Zwart 
|et al.|2010| and references therein). Observational evidence for this 
comes from a number of channels: observations have directly mea- 
sured the global correlations functions of young O-stars (Parker 
|& Goo dwin 2007), and shown that most can be directly identi- 



fied as part of young clusters/associates (Gies|1987 1, with most of 
the remainder being identifiable as runaways ( |de Wit et al.|2005| 
|Schilba ch & Roser 2008 ). The observed star formation rate in em- 
bedded clusters in the MW can account for most of the large-scale 
average SFR ( |Lada & L ada 2003 ). The correlation function of stars 
in various MW star-forming regions has been measured and rises 
continuously on small scales, over > 4 — 5 orders of magnitude in 
radius (Gomez et al.||1993| |Larson||1994[ |Simon||1997l |Nakajima 



et al. 1998 Clarke et al. 2000; Hartmann 2002 Hennekemper et al. 



2008 Kraus & Hillenbrand 2008). In nearby galaxies the obser- 
vations are more difficult, but in starburst galaxies a large fraction 
of all the UV light is often identified with just a few young star 
clusters (see e.g. |Meurer et al.|1995||Zhang et al.|2001| >. 

Theoretically, this is broadly understood as a consequence of 
dense gas being concentrated in primarily in GMCs, which then un- 
dergo fragmentation and turn some fraction of their mass into stars. 
And it does appear the same clustering is evident in the pre/proto- 
stellar core population ( Stanke et al. 2006 1. However, this does not 
provide a quantitative model for their correlation function, nor does 
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it actually explain why this should be the case - why do most stars 
not form in relatively more isolated high-density regions? 

Although much theoretical progress has been made in under- 
standing how stars form in clusters, our fundamental understanding 
of why their formation is clustered remains quite limited. There is 
no analytic theory for the star-star correlation function, outside of 
very small scales where it is dominated by the binary star popu- 
lation. Much of the discussion in the literature has focused on de- 
termining the fractal dimension of star formation: if the ISM were 
structured in a hierarchically fractal manner, this would imply a 
simple power-law correlation function from which the fractal di- 
mension can be determined. But this is not a physical, predictive 
model; it is a parameterization - albeit a useful one - of the obser- 
vations (see also Bate et al. 1998). Moreover it must break down 
at both small scales (the binary regime) and large scales, since it is 
also observationally established that star formation is not strongly 
clustered/biased relative to the dense gas in the disk on very large 
scales ( |Zhang et al.|200T]|Leroy et al.|2008|rFoyle et al.|2010^ 

At first glance, it is not surprising that there is no more so- 
phisticated analytic theory for the correlation function of stars. The 
process is highly non-linear, chaotic, and involves a wide range of 
physics including turbulence, cooling, magnetic fields, and feed- 
back, and the correlation function is a spatially dependent quantity. 
There have therefore been some attempts to compare stellar cor- 
relation functions with numerical simulations of star formation in 
clusters (see Klessen & Burkert 2000; Hansen et al. 2012), but the 
dynamic range of the problem makes this very difficult. Galaxy- 
scale simulations are required to follow the formation of GMCs 
and predict their global properties, and how these overdense re- 
gions compare to the density field in more isolated regions (i.e. 
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why there are not more stars in non-cluster regions). But these can- 
not resolve the formation of stars and model it with sub-grid recipes 
- the typical "star particle" in such simulations is the mass of a real 
star cluster. Simulations which have the necessary resolution can 
only follow small "protocluster" regions whose properties must be 
assumed as initial conditions; they can predict quantities such as 
the binary distribution, but not global clustering. 

However, Hopkins (201 1 1 (hereafter Paper I) recently general- 
ized the excursion set formalism - well known from cosmological 
applications as a means to calculate halo mass functions and clus- 
tering - to calculate the statistics of bound objects in the density 
field of the turbulent ISM. The key property of supersonic turbu- 
lence that makes this possible is that the density distribution (at 
least outside of those bound regions) tends towards a lognormal, 
with a dispersion that varies in a well-defined manner with Mach 



number (see e.g. 


Vazquez-Semadeni|1994||Padoan et al.|1997||Os-| 


|trikeretal.|1999 


i. In Paper I, we use this to construct the statistics 



of the "first-crossing distribution": the statistics of bound objects 
defined on the largest scales on which they are self-gravitating. 
We showed that the predicted mass function and correlation func- 
tions/clustering properties agree well with observations of GMCs 
on galactic scales. In Hopkins (2012) (Paper II), we extended the 
formalism to the "last crossing distribution" - specifically, the mass 
function of bound objects defined on the smallest scales on which 
they remain self-gravitating but do not have self-gravitating sub- 
regions (i.e. are not fragmenting). We argued that these should be 
associated with proto-stellar cores, and in Paper II showed that the 
resulting core MF agrees well with canonical MW CMF and (by 
extrapolation) stellar IMFs. This formalized the approximate argu- 
ment for the same in Hennebelle & Chabrier (20081, but more im- 
portantly for our purposes here, places it in the proper excursion set 
formalism and so makes it possible to calculate higher-order statis- 
tics such as clustering properties, using only information on global 
(galactic) scales. 

In this Letter, we use the excursion set formalism, in anal- 
ogy to its well-studied application for the clustering of dark mat- 
ter halos, to develop a fully analytic a priori theory for the clus- 
tering of protostellar cores and, by extension, star formation, in a 
turbulent ISM. In § [2] we derive the solution for the statistics of 
the last-crossing distribution inside a large-scale over/underdensity 
(the mathematical underpinning of the correlation function). In §[3] 
we show how this relates to the conditional mass function of cores. 
In §[4] we use this to derive the star-star correlation functions. In 
§ [5] we present the calculated correlation functions as a function 
of stellar mass and spatial scale, as well as global turbulent ISM 
properties, and compare to observations. In §[6] we summarize our 
results and discuss their implications. 
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Figure 1. The predicted 3-dimensional auto-correlation function of collaps- 
ing cores, as a function of radius (in units of the disk scale height h). The 
excursion-set model for collapsing cores is determined by solving the two- 
barrier problem for bound/collapsing objects defined at "last crossing" - i.e. 
the smallest scales on which they are self-gravitating. The model is com- 
pletely specified for a given turbulent spectral index p and normalization 
(which we take to be the Mach number at the scale h, Mi,)- Here we take 
p = 2, Mh = 10. Solid lines show the exact solution from Eq. |24| Dotted 
lines compare the approximate scaling in Eq. |28| - the normalization is sys- 
tematically too large, but the shape is correct up to scales near ~ h, where 
the assumption that the run in S(R) is small breaks down. The shape is sim- 
ilar at all stellar masses, but the small-scale amplitude increases strongly at 
low masses. In all cases, cores/stars are predicted to cluster very strongly on 
small scales, below their parent GMC scale. On large scales (> h), they are 
only weakly clustered. The characteristic scale is set by the fact that most 
of the turbulent velocity power (hence power in density fluctuations) is at 
this scale. 
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2 THE TWO-BARRIER LAST-CROSSING PROBLEM 

Calculating clustering properties in the excursion set formalism de- 
pends on the solution two the "two-barrier problem" - the proba- 
bility of an "event" given a specific value of the Gaussian field at 
some larger scale. This is well-studied for the first-crossing events, 
but has not been calculated for the last-crossing distribution; we 
therefore calculate this here. Our derivation closely follows that in 
Paper II. 

Consider the Gaussian field S(x \ R) (which here represents the 
logarithmic density field smoothed in a kernel of radius R about x). 
The variance S(R) is monotonic so we take S as the independent 
variable and consider S(S) (for convenience, we will drop explicit 



Figure 2. The (exact) predicted 2-dimensional correlation function between 
cores of all masses - this should be comparable to the observed (pro- 
jected) star-star correlation function (averaged over the IMF range of stellar 
masses). The predictions are qualitatively similar in all cases, but compli- 
cated second-order differences arise. Models with turbulent spectral index 
p = 5/3 are similar to those with p = 2 but a higher .A/!/, by a factor ~ 2, 
M(r) <x A p ~ ')/ 2 falls off more slowly at the small scales of interest. At 
fixed p, lower values of Mi, produce slightly steeper, more power-law like 
behavior. The more pronounced flattening of £ on small scales with higher 
Mi, arises because the line-of-sight integral is dominated by the (larger) 
power near large scales. 
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Figure 3. Mean surface density E„ (R) of stars around a random star: this 
is proportional, by definition, to 1 + (?2d(S)} in Fig. [2] We compare the 
model predictions (linestyles as Fig. [2] normalized by the mean surface den- 
sity of the observed systems and scale height h = 200 pc) to observations 
of the corresponding core-core cross-correlation (large circles; Stanke et al.| 
(2006) ) and a compilation of observations of the star-star cross-correlation 
(lines and grey points from |Simon] (l 997 1 ; Nakajima et al. (1998i; Hart- 
|mann| (2002) ; |Hennekemper et al.|pOQ8) ; |Kraus & Hillenbrand| (2008) ). 
The agreement is reasonably good over the dynamic range observed, and 
the observational scatter in the shape of §2d is similar to that predicted in 
Fig. [2] The models do not extend to the sharp rise in clustering at scales 
< 0.01 pc, usually attributed to binaries. 



notation of x). The PDF of S(S) is, by definition, 



(1) 



The barrier B(S) is the minimum value S(S) which defines 
objects of interest (e.g. collapsing regions). Consider a "trajectory" 
5(S) which begins at some <5, at 5, (i?, — > 0), then evaluate it at suc- 
cessively larger scales (smaller S). We define last-crossing distri- 
bution in Paper II, fi (S \ 6t) dS, as the probability that the trajectory 
crosses B(S) for the first time between 5 and S + dS without having 
crossed B(S) at any larger 5. The sum over all trajectories is just 
given by f t (S) = / f e (S \ 8,) P(S,) dS,. 

Define fe(S\ So) as the value of fe(S), for trajectories which 
have a value 8 (So) = So at some larger scale So > S. Also define 
n(5[5] I <5o[So]) d<5 as the probability for a trajectory - constrained 
to have 5o(So) - to have a value between 8 and 8 + dS at scale 5, 
without having crossed B(S) at any larger S' > 5. The integrals of 
fe(S I So) and n(5[5] | (JoISo]) for 5 < B(S) must sum to unity: 



KS) 



dS'f e {S'\S )+ Il(6[S]\So[So])d8 (2) 



1 = 



If we ignored the barrier, II(<5[S] j <5o[5o]) would be equal to 
/'io(5[S] |<5o[5o]), the probability that a trajectory with 8o(So) has 
a value 8(S) at a larger S > So- This is just 



Pw(6[S] I 80IS0]) = P (S-S \S- So) 



(3) 



But we must subtract from this the probability that a trajectory 
crosses the barrier at some larger S' > S and then passes through 



5(S), so 

n(5[5]|5o[5o])=Pio(<5[5]|5o[5o])- 



f 

■Is 



dS' f e {S'\5,)Poi[5(S) \S' =B(S'), So(So)] 



(4) 



Where Poi [8(S) | S'(S'), So(So)] is the probability of a transition 
from 8(S') to 8(S) where S' > S - i.e. moving in the "opposite" 
direction (decreasing 5) from the transition which defined Pw - for 
a trajectory subject to the constraint that it must equal 80 at So- This 
can be determined from Bayes's theorem: 

P l0 (SlS]\8 [S ]) 



Poi[8(S)\B(S')]=P lo [B{S')\8(S)]- 



(5) 



P 10 (S[5']|<5o[5o]) 

where on the right-hand side we have Pio(S[S] | <5o[So]) and 
Pio(B[S'] I 5 [So]) instead of P (8 1 5) and P (B[S'] j 5") because the 
distributions are constrained to have <5o(5o). But because fluc- 
tuations on successive scales are uncorrelated, the probability 
PiolB(S') j S(S)] does not explicitly depend on So- 

The governing equation for fe(S\ So) is now completely de- 
fined in terms of normal distributions. We can therefore follow ex- 
actly the procedure in Paper II. After differentiating Eq.|2]to isolate 
fi (S) and performing a considerable amount of simplifying integral 
evaluation, we arrive at the key equation for ft (S \ So) : 



MS\So)= gi (S\S )+ dS' fe (5' \S ) g2 (S, S'\S ) (6) 



where 



gi(S\So) 
gi(S, S'\S ) 



, dB B{S) - So 



dS S-So 
B{S')-B(S) B(S)~Sp 
S'-S S-So 

P [B(S)-5 e (S, S')\S e {S, S') 



]P (B(S)-S \S-S ) (7) 
(8) 



- 2 ds\ x 



5e (S,S')=B(S') + (S -B[S'l)(^-£j 
S e (S,S') = (S'-S)( s ~ s " 



(9) 
(10) 



^S'-So 

Although complicated, Eq.[6]has, in general, a unique solution for 
any arbitrary barrier B(S), and is straightforward to solve via stan- 
dard numerical methods for any choice 8o(So) provided So < S (see 
the discussion in Paper II and Zhang & Hui (2006)). 

In fact for a linear barrier, B(S) = Bo + /3 S, this has a remark- 
ably simple closed-form solution: 

fe (S I 80) \b=b +ps = PPo{B[S - So] -5o\S- So) (11) 

It is easy to verify that when So — > as So — > 0, we recover 
the one-barrier last-crossing distribution from Paper II, which is the 
value of fe (M j So ) averaged over all So (So ) 

fe(M) = (f e {M I So)) = fe(M S -> 0, So -> 0) (12) 

3 THE CONDITIONAL MASS FUNCTION 

In Paper I we derive S(R) and B(S) from simple theoretical consid- 
erations for all scales in a galactic disk. For a given turbulent power 
spectrum, with the assumption that the disk is stable (Q — 1), S(R) 
is determined by summing the contribution from the velocity vari- 
ance on all scales R' > R 



S ( R )=f o |W(/:,^)| 2 ln[l + 



4 c* + K 2 k- 



dlnk 



(13) 
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where W is the window for the density smoothing[^|s(i?) is 

' Pcrit ^ , S(R) 



B(R) = In 



Po 



Peril 
PO 



(14) 



(15) 



where po is the mean midplane density of the disk, ft = k/Q, = v2 
for a constant- V t . disk, and 



oJ(*) = c? + vi + (v?(U)) 
The mapping between radius and mass is 



M(R) : 



r R 

^P^[— 2 + 



l + : 



R 



exp 



(16) 



(17) 



It is easy to see that on small scales, these scalings reduce to the 
Jeans criterion for a combination of thermal (c s ), turbulent (v,), and 
magnetic (va) support, with M = (47r/3) pcniR 3 ', on large scales it 
becomes the Toomre criterion with M = irYi a nR~ . 
The mass function is given by 



dn 
dM 



Pc(M) 
M 



MM) 



dS I 
dwl 



(18) 



Two parameters completely specify the model in dimension- 
less units. These are the spectral index p of the turbulent veloc- 
ity spectrum, E(k) oc k~ p (usually p « 5/3 — 2), and its normal- 
ization, which we define by the Mach number on large scales 
M\ = {v 2 (h)} I (c 2 + Va), The dimensional parameters h (or c s ) and 
po simply rescale the predictions to absolute units. 

For a choice of p and Mi,, it is straightforward to numeri- 
cally determine ft (M | <5o[.Ko]), hence the conditional mass function 
(mass function within the appropriate sub-regions for any choice 
of So and Ro). Unfortunately a closed-form solution is not gener- 
ally possible. However, we show in Paper II that on sufficiently 



small scales (near/below the sonic length W SO nic = hM 



-2/0,-1) 



). 



the "run" in S(R) ~ S(R SO mc) becomes small (since most of the 
power contributing in Eq.[T3]comes from large scales), while B(R) 
rises rapidly, so dB/dS 3> B(S) /S^> 1 . In this limit, the conditional 
MF is approximately 



dn 
dM 



Peril 

(M) Idlnp 

CI 



M 2 



dlnM 



\Po(S-Sa\S-Sa) 



where 



Po(S-S \S-S ) 



1 



exp 



V2tt(S-5o) 
(ln[p cnt /p(5 )] + (5- 



(19) 



(20) 



So)/2) 2 



2(S-S ) 

For the appropriate choice of So on scale Ro, this is the result- 
ing core MF (and so relates to the resulting stellar IMF and SFR) in 
over-dense regions (e.g. GMCs) or "voids" (inter-GMC gas). Tak- 
ing So — > and Ro — > oo, we recover the galaxy-averaged CMF. 

4 THE CORRELATION FUNCTION 

The auto-correlation function £ of a given population is denned as 
the excess probability of finding another member of the population 

' For convenience we take this to be a &-space tophat inside k < 1 /R, which 
is implicit in our previous derivation, but we show in Paper I and Paper II 
that this has little effect on our results. 



within a differential volume at a radius r from one such member, 
i.e. 



(N{r\M)) 
(n(M))dV 



(21) 



But since n(r\M) oc fe(M\ 8o[So{r)\), for a region with an over- 
density So on a scale So (corresponding to r), this is just 

l + fMM(r|M)=/ ( M f M l S . o) )p t (So\So[r])d8 (22) 
J s v Je( M ) ' 

Here, P, (So \ So[r]) is the probability of 5 (So) having the value So 
on the scale So, given that 8(S[M]) = B(S[M]) ~ i.e. that there is a 
barrier crossing (a core/star) at the "starting point." But the proba- 
bility of having a last-crossing event at scale S > So, given a density 
<5o(So), is just fe(M \ So). So conversely, the probability of <5o(So) 
given that crossing is just related by Bayes's theorem 

^l^ ^f'^^ /^l^^g (23) 



So we obtain 



(fla S .cro ss (M)) ' f e (M) 

l+Zuu(r\M) = f°° ( M f M l 5 ° ) ) 2 Po(So\So[r])dSo 

i-oo v M M ) ' 



(24) 



More accurately, what is typically measured is the star-star 
or clump-clump correlation function over a broad mass range. We 
therefore require the cross-correlation between an initial crossing at 
Mi and second crossing at some Mt. The logic is identical, however. 
We then integrate over the appropriate range of Mi and, finally, 
average over Mi (weighted by number density). We obtain 



i+^W 



where 



/ dW Ml J dW Ml ft (Mi \So)fi {Mi | So) P(S | So) dS 



(J dW M fe(M)) 2 



dW M 



Pent (M) |dS(M) 



M 



dM 



dM 



(25) 



(26) 



and the integrals over mass M should be over the appropriate ob- 
served core/stellar mass range. 

Finally, note that what is generally measured is the projected 
correlation function £id(R p ). Projecting £ (r) is straightforward: 



tufa) 



JZo"°(z)dz 



(27) 



where z is the line-of-sight direction and no(z) is the average abun- 
dance expected. Technically this will depend on projection angles 
(especially in the complicated case of systems inside the MW), but 
for almost any profile where no ~ n(M) is a weak function of z and 
then falls off at scales > h the results are nearly identical. 

5 RESULTS 

In Fig. [T| we plot the exact (3d) autocorrelation function from 
Eq. [24] as a function of r and M. Recall, in dimensionless units 
the problem is completely specified by a choice of p and Mh, for 
which we adopt canonical values p = 2 (Burgers turbulence, typ- 
ical in the supersonic regime) and Mi, — 10 (typical of the cold 
clumps in the MW disk). 

The exact results in Fig. [T] require the numerical solution for 
fi(M \So). However, if we restrict to the regime where dB/dS 2> 
B/S 3> 1 so that the mass function can be approximated by Eq. 19 
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or consider a linear barrier (Eg. 1 1 1), then this has the closed-form 
solution 



l+£mt(r\M) ! 



: exp 



v/l-(So/S) 2 ~' r \S(l+S/So) 
At small and large r, respectively, this becomes 
exp(B 2 /2S) 



(mm('" I M) 



£,uu{r\M) 



V2(l-5 /S) 
(Sb- 



(So -» 5) 



0) 



(28) 



(29) 



(30) 



We compare this to the exact result. The shape of £(R) and its sys- 
tematic dependence on R are qualitatively described. However, it 
is nowhere exact, because S(R) does vary significantly on larger 
scales, and B(R) and S(R) both vary such that even if dfi/dS S> 1, 
dB/dS > B/S is less true. 

The shape of £mm(»0 is similar in all cases - i.e. it is nearly 
mass-independent over the interesting range - so even though the 
expression for the CMF-averaged cross-correlation £(r) in Eq. 25 
is complicated, it mostly amounts to a MF-weighted normaliza- 
tion. What we care about is the correlation function shape (we note 
below that the observed normalization is arbitrary and should be 
factored out in any case). This means that regardless of the CMF 
shape used to "average" the prediction, the result is similar. For the 
same reason, provided the general assumption that stars form from 
self-gravitating cores is reasonable, there will be no large differ- 
ences between the predicted stellar and core correlation functions, 
regardless of how this process occurs. If the conversion from core 
mass to stellar mass were constant, the correlation functions would 
be exactly identical; but even allowing for a large random variation 
in conversion efficiency or systematic mass dependence makes lit- 
tle or no difference to our conclusions. We have, for example, con- 
sidered sampling Eq.[25]over only different sub-ranges in mass (say 
sampling only massive stars), or sampling over the full stellar IMF 
but with a random 0.5 dex scatter in the assumed stellar-to-core 
mass ratio, or allowing for the stellar-to-core mass ratio to depend 
on mass oc M" 1 ; in all cases the differences in the predicted f (r) 
shape are smaller than the differences owing to different choices 
of the turbulent spectrum. What will change the clustering, on very 
small scales <C R i0 mc, is if individual single Jeans-mass cores form 
multiple stars. And to some extent, this must happen, as many stars 
are in binary systems. This is not accounted for here; therefore the 
predictions cannot be reliably extended to scales below the charac- 
teristic core scale < R SO nic (and indeed we see a discrepancy appear 
at these scales). We discuss this further below. 

In analogy to the clustering of dark matter halos and GMCs 
defined in Paper I, we can define the linear bias as the ratio of the 
autocorrelation to the autocorrelation function of the mass itself. 
On large scales, this should approach a constant, 



b(M) 2 = £/^ : 



('■- 



(31) 



On large scales (where So[r] is small), the autocorrelation of the 
mass is just the variance £ mass « So (see |Mo & W hite 19961. From 
Eq. [30] we can therefore immediately obtain the bias in the limit 
dB/dS>B/S>l, 



b(M) : 



B(S) 
S 



1- 



0.04 In (M/M S( 
0.17 In (M/M s< 



M > M s< 



(32) 



where the latter equality is for typical parameters (p = 2, Mh ~ 30) 
and Msoniq = M(R son i t: ). The large-scale bias is a very weak function 
of mass and about unity. 

In Fig. [2] we plot the projected ^2d from Eq. [27] for the full 



CMF/IMF-averaged £(r)[]We examine the effects of the two free 
model parameters: the large-scale Mach number Mi, and turbulent 
spectral index p. As shown in Paper II for the CMF/IMF, chang- 
ing p — 5/3 instead is nearly equivalent to assuming a higher Mi, 
at p = 2 (since, in both cases, the sonic length is pushed down and 
turnover at low masses is slower). At lower Mi,, decreasing Mi, in- 
creases the clustering signal, because there is a smaller variance 5, 
so the small-scale cores depend more heavily on large-scale fluc- 
tuations. At sufficiently large Mi,, however, there is a significant 
change in shape in ^2d; the initial rise below ~ h is rapid, but near 
^sonic, £.26 flattens. This is a projection effect: the rise in £34 on small 
scales is still similar (just slightly more shallow) to that we see in 
£2d at smaller Mi,\ however there is much more power on inter- 
mediate scales as well, so the power in ^2d is nearly "converged," 
hence the flattening. 

In Fig. [3] we compare the projected £ia to a compilation of ob- 
servations of both the star-star and core-core autocorrelation func- 
tions. Typically, the observations are plotted not as ^2d but as E„, 
defined as the average surface number density dN/dA(i? p ) of com- 
panions at a radius R p from the primary. By definition, this is triv- 
ially related to ^2d as 



E.(*,) = (E.)(l + 6d) 



(33) 



Here, E* is undetermined - it depends both on the absolute back- 
ground galaxy properties which set the scale of the problem (h, 
po, etc) and on the uncertain star formation efficiencies (fraction of 
gas which will actually turn into stars). We treat it as arbitrary and 
simply compare the profile shape. We normalize R to an absolute 
physical scale by assuming a MW-like scale height h = 200 pc for 
the gas. 

We can gain some (approximate) insight into the functional 
form of ^2d on small scales from Eq. [28] Assume r <C /i (the dy- 
namic range of most interest), in the limit of Eq. 



29 



so M\r) < 

Mf, as well, and expand to leading order in O(rfh) and 0(M^ ) 
(series expanding Eq. |13| around 5 = So); finally make the simple 
approximation that projection to £>d multiplies £(r) by one power 
of R and that we can treat the mass-averaging as a normalization 
correction. After some tedious algebra, we obtain 



1 + f 2d OC 



(l + |^-')'/4 + (^) (1 +x „-l , 

x("-0/2-l L y/2 



X 2 J 



(34) 



cc (1+x)(1 + ^ )1/4 +0(ln-'M,,) (p = 2) (35) 

where x = R/R sonic , and e = ln(jr 2 (1 +x p ~ 1 ))/4hi(M h ). Al- 
though we have made a number of simplifications, we recover 
all the key behaviors in Fig. [2] the predicted ^2d(R) is an ap- 
proximate power-law (because p a it(r) is so), which scales steeply 
(oc R~^- 5 ~ 2 - 3 ' for p = 5/3 — 2) on small scales, then flattens around 
the sonic length (first to oc r-( - 5 - 10 ) then R~ 025 ). The "steepen- 
ing" below i? son i c is caused by the steeper dependence of p cr i t on R 
at these scales (thermal support preventing collapse); the gradual 
run at larger scales from the logarithmic run in S with M(r). 

Although it is less precisely defined, we can also make some 



2 Technically, we integrate the CMF range from 10" 4 - 10M somc , which 
for typical parameters corresponds to ~ 10 — 3 — 100 Mq , but we stress that 
so long as we include the "peak" of the CMF/IMF, this choice makes lit- 
tle difference. To project, we also assume the average abundance «o(z) 
is flat out to ±2 h, but changing this limit or assuming an exponential 
no oc exp(— z/h) has only weak effects. 
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estimate of the fraction of stars formed in an "isolated" (non- 
cluster) mode. Specifically, in analogy to our derivation of Eq.|24| 
we begin with the probability that a core (last-crossing event) at 
fe(S[M]) is embedded in a larger region So[r] with some <5o(So); we 
can then calculate the probability that this region contains zero ad- 
ditional last-crossing events. This latter probability is just given by 

i.e. one minus the total probability 
of a crossing at 5 > So- After some algebra, we obtain the probabil- 
ity of a core of mass M having no additional crossings within the 
parent radius r[S]: 



M<r\M) = l- 



dS 



: ds 'f^j^l M s'\S )P (So\S [ r ]) 
HV>) 

(36) 



3 IdS I 
-0.15 [in 



If we take the limits dB/dS 3> B/S 3> 1, we can obtain the ap- 
proximate scaling / iso (r < n\M[r ]) ~ \n~\n/r ) (n/r )~ d ' !/d5 . 
In greater detail, if we solve this for a linear barrier and and assume 
the less restrictive dB/dS > 1 and r <^h, and insert B(r) and S(r) 
defined above, we finally obtain 

/iso(<r|5)w||^|~ I Po(B[51-B[5o]|S-5 [r]) (37) 

r M-'/2/ r N-0-65 

D J I (n— ) (38) 

"sonic ' J v J^Nonic ' 

where in the latter we insert p = 2 and Mi, ~ 5 — 30 (the result is 
very weakly dependent on Mi,). This should translate to the frac- 
tion of cores formed without a companion core inside of r. We cau- 
tion that if only a small fraction of cores make stars, the fraction of 
stars formed without another star inside < r will be higher. But if 
cores produce multiple stars (recall that we do not explicitly treat 
binaries here), the fraction of isolated stars will be even lower. Re- 
gardless of these uncertainties, we see that the absolute number 
of cores/stars formed without neighbors inside a radius r is small 
even for r near the sonic length, and falls rapidly as we expand the 
"search radius." 

6 DISCUSSION 

We have developed an analytic theory for the clustering properties 
of protostellar cores in a supersonically turbulent density field. In 
Paper I, we developed an excursion-set theory for lognormal den- 
sity fluctuations in the ISM, and applied it to the mass functions and 
clustering of GMCs - the "first-crossing distribution" (distribution 
of bound masses on the largest self-gravitating scales). In Paper II, 
we showed that this could be extended to predict the protostellar 
core MF by defining the "last-crossing distribution" (distribution 
of masses on the smallest scales on which they are self-gravitating 
and non-fragmenting). Here, we extend this method to define the 
conditional last-crossing mass function, as a function of density on 
larger scales: i.e. the "two-barrier last crossing problem." We use 
this to analytically derive the correlation function of these cores as 
a function of separation (and core mass), on all scales from the core 
itself to galactic scales. This should directly translate to the corre- 
lation properties of newly-formed stars, independent of the mass 
conversion/star formation efficiency. 

We find that, on large scales > h (the disk scale height), cores 
are weakly biased. Star formation, therefore, should not be es- 
pecially sensitive to galaxy-scale overdensities, except insofar as 
those overdensities collect/drive the dense gas. This has been ob- 
served: star formation does not appear to be strongly biased rela- 
tive to the gas itself in spiral arms, etc. (Leroy et al.|2008]|Foyle| 
|et al.|2010) . The reason is simple: star formation is a primarily lo- 
cal process. It is also closely related to our derivation in Paper I of 



the clustering of GMCs: they themselves are not strongly biased 
(relative to the rest of the dense gas) on large scales, provided there 
is sufficient gas present and cooling is efficient. 

On scales <C h, however, cores cluster very strongly. The for- 
mal three-dimensional clustering amplitudes at small scales are 
enormous - in other words, stars form in clusters. The approximate 
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illustrates 



13 (hence 



conditional mass function (on small scales) in Eq. 
why. On small scales, the "run" in S(R) given by Eq. 
A5 = 5 — 5o in Eq.|l9[ is small, because M 2 (R) is small. As this 
vanishes, the Pq term in Eq.[l9]becomes extremely sharply peaked 
around p(Ro) = Pmt(R) - so the CMF/IMF-averaged conditional 
mass function approaches a step function in p(Ro) on some larger 
scale. Given the same arguments, the fraction of stars formed in an 
"isolated" mode should scale as O (dB/dS) ~ 1 ; near the sonic length 
this is approximately \dB/dS\^l R < 0.1. 

Physically, the meaning of this is: as we average on smaller 
scales approaching the sonic length, the local rms Mach numbers 



decline rapidly (M 



"'), so the contribution to density fluctua- 



tions (the variance S(R), which is directly tied to the Mach number 
in supersonic turbulence) also drops. But the threshold density for 
collapse continues to rise. So the smallest scale at which a region is 
going to be self-gravitating (or avoid fragmentation) is essentially 
"imprinted" by the density by fluctuations driven on larger scales. 
It is unlikely to "wander" much across the barrier on smaller scales 
since fluctuations are small. 

As a result, stars preferentially form in larger-scale overdense 
regions which must themselves (on some scale) be self-gravitating. 
The characteristic scale of the "parent" systems - i.e. the charac- 
teristic scale at which the clustering amplitude will fall off, is the 
scale where S(R) begins to run significantly with R. But this is by 
definition the characteristic scale of the first-crossing distribution, 
which we identified in Paper I with GMCs, and showed there has a 
characteristic scale of ~ h — just set by the global turbulent leans 
length. In short, stars form inside of GMCs, in strongly clustered 
fashion, with the characteristic length of that clustering set by ~ h. 

These arguments are quite general: the fundamental statement 
is that, since most of the power in the turbulent velocity field is on 
large scales (~ h), star formation must be strongly clustered below 
that scale. This is true for any plausible turbulent power spectrum. 
We show explicitly that a qualitatively similar scaling results for 
turbulent spectral indices p ~ 5/3 — 2, and large-scale Mach num- 
bers Mi, ~ 5 — 30. However the details of the correlation function 
do depend on these values, in a non-trivial manner. For example, the 
slope on small scales ^2d oc B~ a , is an approximate power-law with 
values a ft* 0.5 — 1 .0; it can flatten significantly to a ~ 0. 1 — 0.5 on 
sufficiently small scales if the Mach numbers are sufficiently large. 

We compare the predicted correlation functions to observa- 
tions of both protostellar cores and young stars, and find that they 
agree well on the applicable scales. The scatter in the shape of ob- 
served correlation functions is also intriguingly similar to the range 
predicted from the parameter variations above. 

On the smallest scales <0.1 pc, our prediction does not rise 
as steeply as the observed stellar correlation functions. This is ex- 
pected, since this effect comes primarily from binary stars, which 
we do not explicitly predict since the fragmentation events are on a 
sub-core scale. It may however, be possible to treat this in the fully 
time-dependent formulation of the excursion set ISM model devel- 
oped in Paper I, since that could follow the simultaneous growth of 
fluctuations and collapse of a given sub-region. 
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